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Abstract 

We focused our attention on Ca 2+ release from the endoplasmic reticulum through a clus- 
ter of inositol 1,4,5-trisphosphate (IP3) receptor channels. The random opening and closing 
of these receptors introduce stochastic effects that have been observed experimentally. Here, 
we present a stochastic version of Othmer-Tang model for IP 3 receptor clusters. We address 
the average behavior of the channels in response to IP 3 stimuli. We found, by stochastic 
simulation, that the shape of the receptor response to IP3 (fraction of open channels versus 
[IP3]), is affected by the cytosolic Ca 2+ level. We also study several aspects of the stochastic 
properties of Ca2+ release and we compare with experimental observations. 
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1 Introduction 



Ionized calcium (Ca 2+ ) not only represents the most common signal transduction element relay- 
ing information within cells to control a wide array of activities including secretion, contraction 
and cell proliferation, but also is invariably involved in cell death (Berridge, 1997). To coordinate 
all of these functions, cytosolic Ca 2+ needs to be precisely regulated in space, time and ampli- 
tude. Normal concentrations of cytosolic Ca 2+ ([Ca 2+ ] ~ 100 nM) is 20000 fold lower than the 
2 mM concentration found extracellularly. Under resting conditions this gradient is maintained 
by active extrusion of cytosolic Ca 2+ by Ca 2+ ATPases present in the plasma membrane and in 
the endoplasmic reticulum (ER), or in the sarcoplasmic reticulum in electrically excitable cells. 
These pumps counterbalance the leak of Ca 2+ into the cytosol from both the extracellular space 
and the ER Ca 2+ store. 

A wide variety of extracellular stimuli cause the increase of [Ca 2+ ] to exert their effect. In 
non-excitable cells, this increase is triggered by inositol(l,4,5)-trisphosphate (IP3), produced 
upon activation of phospholipase C (Berridge, 1997). IP 3 rapidly diffuses into the cytosol, 
where it interacts with the inositol(l,4,5)-triphosphate receptors (IP3R), to promote the release 
of Ca 2+ into the cytosol. Depending on the cell type, the resulting cytosolic Ca 2+ signal can 
have a complex spatio-temporal composition (Berridge, 1997). It is generally recognized that, 
in the presence of a constant external stimulus, the Ca 2+ displays spiking behavior. These Ca 2+ 
oscillations can be spatially localized or extended, spreading as a wave throughout the entire 
cell. 

The process of Ca 2+ released from the ER through channels is nonlinear since, increased 
Ca 2+ concentration in the cytosol favors channel opening. This autocatalytic amplification is 
called calcium-induced calcium release (CICR). There are a variety of channels showing CICR. 
Ca 2+ release is terminated by closure of the channels at high Ca 2+ levels, after which Ca 2+ is 
removed from the cytosol by the action of the Ca 2+ ATPases. 

It has been observed that Ca 2+ release channels are spatially organized in clusters. Single 
release channel, named Ca 2+ blips, has been observed in experiments (Bootman et al., 1997). 
Collective opening and closing of several Ca 2+ channels in a cluster, named puffs, have also 
been observed in experiments (Callamaras et al., 1998). That suggests a hierarchy of calcium 
signaling events from small blips to large puffs (Lipp and Niggli, 1998; Bootman et al., 1997). 
Improved spatial and temporal resolution in recording reveal that there is not a clear distinction 
between fundamental blips and elementary puffs. Channels open and close in a stochastic way. 
Ca 2+ release by one channel increase the open probability for the neighboring channels. Thus 
at high levels of IP3, neighboring clusters become functionally coupled by Ca 2+ diffusion and 
CICR support the formation of spatiotemporal patterns of intracellular Ca 2+ release. 

Watras et al. (1991) found that Ca 2+ channels in the cerebellum exhibits four conductance 
levels that are multiples of a unit conductance step. This observation suggests that the number 
of interacting receptors in one complex can vary from one to four and supports the hypothesis 
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that the channel is a tetramer. Examination of the IP3 dependence of channel, yielded Hill 
coefficients of 1-1.3 (Watras et al., 1991). However, in hepatocytes the response to IP3 are 
positively cooperative with a Hill coefficient of 3.0-3.4 (Marchant and Taylor, 1997; Dufour et 
al., 1997) or. In our numerical simulation we found that the Hill coefficient increase from 1.15 
to 2.2 when total Ca 2+ (cytosolic + stored) level increase. However, the model is not able 
to explain higher values of as that reported by Marchant and Taylor (1997). This finding 
suggest that feedback from cytosolic Ca 2+ plays a key role in the channel response to IP3, but 
is not able to explain high Hill coefficients found in some tissues. 

Experimental finding of Bezprozvanny et al. (1991) clear up several aspects of the kinetic 
of the receptors. Thus several models for Ca 2+ dynamics in systems involving IP3R have been 
proposed (De Young and Keizer, 1992; Othmer and Tang, 1993; Bezprozvanny and Ehrlich, 
1994). These models of IP3R assume a regulatory site for IP3 on the channel complex, one 
activating site and one inhibiting site for Ca 2+ . Experimental findings suggest that the channel 
is open when IP3 ion is bound to its corresponding domain and Ca 2+ is bound to its activating 
domain and not bound to its inhibiting site. In the aforementioned models, Ca 2+ binding to 
the activating regulatory site is a fast process compared with that of inhibiting binding and all 
these models reproduce a bell-shaped curve when one plots the fraction of open channels as a 
function of log[Ca 2+ ]. 

De Young and Keizer model (DKM) for IP3R assumes that the ligands can bind to any 
unoccupied site on the IP3R independently of the binding status of the other sites, which leads 
to eight possible states of the receptor (De Young and Keizer, 1992). A more general scheme is 
to assume state-dependent binding to each site. In this framework, the simplest model was de- 
veloped by Othmer and Tang (1993) which assumes that the binding order of ions and molecules 
to the receptor is not free. In the Othmer- Tang model (OTM) it is assumed that the binding 
process is sequential (Othmer and Tang, 1993). Ca 2+ binds to the activating site on the receptor 
only after IP3 has bound, and that the binding of calcium to inhibitory site occurs only after 
calcium is bound to the activating site. This sequential binding leads to four possible states 
for the receptor. Bezprozvanny and Ehrlich (1994) proposed a variation to the OTM which 
incorporates an additional closed state corresponding to activating site occupied which decay 
rapidly to the open state. 

The bindings of IP3 and Ca 2+ to the regulatory sites are stochastic events rendering the 
opening and closing of the channel a stochastic process. The small number of calcium channels 
in a cluster indicates that deterministic models might be insufficient. In fact, abortive waves 
cannot be understood in terms of deterministic models, since in these models an excitation 
travels steadily if it travels at all. The stochastic effects are relevant for modeling Ca 2+ wave 
propagation. Furthermore, the observation of localized stochastic Ca 2+ puffs and the rather 
small number of channels creating the localized event, suggest that it is mandatory take into 
account the binding processes as stochastic events. The stochastic dynamics of clustered IP3R 
has been studied by Swillens et al. (1999) which model has 14 states. There are also several 
studies considering stochastic dynamics focusing on the onset of the saltatory propagation of 
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calcium waves due to intercluster diffusion of Ca 2+ (Keizer et al., 1998; Falcke et al., 2000). Shuai 
and Jung (2002b) have been study the statistical properties of Ca 2+ release in a stochastic and 
simplified version of the DKM. 

In this paper we regard a Markov-stochastic version of the OTM to study the stochastic 
properties of the calcium release of IP3R clusters. We decided to stick to the Othmer and Tang 
approach because the OTM has the least number of channel states and kinetic parameters, yet 
it adequately explains the experimentally observed channel kinetics. After a brief description 
of the model in Section 2, we discuss the role of cytosolic calcium and IP3 in regulating the 
fraction of open channels. We show that the fraction of open channels versus IP3 concentration 
fit a Hill curve which coefficient increase with the average intracellular calcium concentration. 
The amplitude, inter-puff interval (IPI) and size distribution of calcium puffs are also discussed. 
Finally, we show that there are no correlation between two consecutive IPIs, but there are some 
correlation between the puff duration and the puff amplitude. 



2 Materials and Methods 
2.1 Stochastic version of OTM 

In the Othmer and Tang model there is only one subunit/receptor in a each channel. The three 
binding sites on the receptor and the sequential binding scheme allows four different states for 
the receptor Xi,j,k- The index i stands for the IP3 site, j for the activating Ca 2+ site, and k for 
the inhibiting Ca 2+ . An index is 1 if an ion is bound and if not. The state transition scheme 
is given by the pathway 

fci[IP 3 ] fc 2 [Ca 2 +] fc 3 [Ca 2 +] 

Xo,o,o — Xi i0 ,o =^ Xi 5 i i0 " Xi,i,i> 

k-i /c_2 fe_ 3 

in which the ki, i = ±1, ±2, ±3, are the rate constant of the state transitions. The channel is 
only conductive when the receptor is in the state Xi^o- The channels are assumed to be close 
enough so that Ca 2+ concentrations can be considered homogeneous throughout the cluster. We 
neglect Ca 2+ diffusion between cluster and environment without accounting for spatial aspects 
of formation and collapse of Ca 2+ waves. We will not consider calcium transport across the 
plasma membrane. Thus the Ca 2+ simply flow or is pumped back and forth between the ER 
and the cytoplasm, and the total intracellular calcium concentration is constant. Thus, the 
dynamic of Ca 2+ in the cytoplasm is governed by the equation 

d [Ca 2+ ] _ 

^ — ^channel ~\~ Jleak J pump t (1) 

where J c hannei is the calcium flux from the ER to the cytoplasm through the IP3R channel, 
Jpump is the calcium flux pumped from the cytoplasm into ER, and Ji ea k is the calcium flux 
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being leakage from the ER to the cytoplasm. The expression for the fluxes are given by 

J channel = V r Jl ( [Ca 2+ ] ER - [Ca 2+ ] ) 

Jleak = « r 7o([Ca 2+ ] ER - [Ca 2+ ]) 
[Ca 2+ ] 4 

J„ - Pi [Ca2+] 4 + ^. ( 2 ) 

where v r is the ratio of the ER volume to the cytoplasm volume, 71 is the maximal Ca 2+ fluxes. 
Nopen is a random variable that represent the number of channels open (receptors in the state 
Xi^o) and N is the total number of channels in the receptor. 70 is the basal permeability of 
the ER membrane in the absence of IP3. For the study of Ca 2+ dynamics in many cell types, 
Ca 2+ exchange with the extracellular medium is much smaller than the Ca 2+ flux across the 
ER membrane (Koster et al., 1993; Thastrup, 1990). For this reason, Eqs. (1-2) can be sim- 
plified by using the volume average intracellular calcium concentration [Ca 2+ ] ai)e =([Ca 2+ ]+v r 
[Ca 2+ ]ER.)/(1 + v r ). The [Ca 2+ ] ai)e is a control parameter whose value we can choose, but that 
it is not a dynamical variable. We then define C =[Ca 2+ ]/[Ca 2+ ] a „ e and rewrite the Eqs. (1-2) 
in the form 

C A 

C = {1 + v r ) ( 7o + 71 N open ) (1 - C) - Pi 4 , (3) 

u + p 2 

where p\ = Pi/[C& 2+ ] ave and p 2 = P2/[Ca 2+ ] a „ e . 

This equation was numerically integrated using fourth order Runge-Kutta algorithm, with 
the actual configuration of open channels N open , which is obtained from stochastic simulations 
of the kinetics of the receptors. The number of open channels changes due to six stochastic pro- 
cesses: the IP3 binding (unbinding) to (from) receptors; the activation (deactivation) of receptors 
by binding (unbinding) Ca 2+ to (from) their activation domains; the inhibition (activation) of 
receptors by binding (unbinding) Ca 2+ to (from) their inhibition domains. Transitions involving 
binding (but not those involving unbinding) depend on the concentration of IP3 or Ca 2+ . Let 
us define the probability distribution vector V (t) = (-Po,o,o {t) , A,o,0 (t) > Pi, 1,0 (t) > ^1,1,1 (£)) T > 
where Pij,k denotes the probability that the receptor is in the state Xjj^. The dynamic of the 
vector V (t) is described by the master equation, 

V = QV, (4) 
where Q represents a (4 x 4) matrix of transition probabilities, 



Q = 



(5) 



/ 1 - h [IP 3 ] fc_i \ 

h[lP 3 ] 1 - - A; 2 [Ca 2+ ] fc_ 2 

fc 2 [Ca 2+ ] l-^ 2 -fc 3 [Ca 2+ ] fc_ 3 

V fc 3 [Ca 2 +] l-fe_ 3 y 

For the stochastic part of the simulation, the state of each receptors at time t + At were in- 
dependently updated from the previous state and Q. If the time step At is sufficiently small 
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to allow at most one transition, the probabilities for the receptor to transform into the various 
possible states within this short time interval are given by 

Po,o,o->i,o,o = ki [IP 3 ] At Pi,o,o-o,o,o = k^At 
Pi,o,o-i,i,o = k 2 [Ca 2+ ]At Pi,i,o-i,o,o = k -2 At 

^1,1,0-1,1,1 = k 3 [Ca 2+ }At ^1,1,1-1,1,0 = fc-sAt. (6) 

For this model, there is only one possible transition out of the states Xo,o,o and Xi^i, while 
there are two possible transitions out of the states Xi ; o,o and Xi^o. In the last case a random 
number p between and 1 was drawn from a uniform distribution for each update step. Then, a 
channel initially in state, for example, Xi^o, will be updated to state Xi 1 0, if p < Pi 0,0— >i,i,0j 
and to state X 0)0 ,o if Pi,o,o->i,i,o < P < Pi,o,o->o,o,o + Pi,o,o-i,i,o, and remained in its current 
state otherwise. 

All over the paper we have used A = 0.001 s, and the number of receptors in a cluster, 
N, was setted to 20. This number is based on the theoretical predictions of the cluster size 
by Swillens et al. (1999), from numerical studies regarding the requirements of interchannel 
communication. Furthermore, found that the optimal cluster size for coherence resonance is 
around N = 20 (Shuai and Jung, 2002a). The values of the parameters for the kinetic model 
used in this paper are given in Table 1. 

To observe puffs in experiments, calcium diffusion is suppressed by intracellular loading 
with the Ca 2+ buffer EGTA (Thomas et al., 1998; Marchant et al., 1999; Cheng et al., 1999; 
Callamaras and Parker, 2000). With a large loading of EGTA, the clusters become functionally 
isolated (Callamaras and Parker, 2000). Under these condition, our model is valid. 



3 Results 

3.1 The channel response to IP 3 and Ca 2+ 

Before to study the statistical properties of the intracellular Ca 2+ release, we firstly address 
the average behavior of channels in response to IP3 stimuli. In this sense, we simulate the 
stochastic dynamics (Eq. 4-5) with the parameters in Table 1 and using [Ca 2+ ] am as control 
parameter. We compute the mean value of cytosolic Ca 2+ concentration in a long time period 
(10,000 seconds) as a function of [IP3] for several values of [Ca 2+ ] sm (0.4, 0.8, 1.2, 1.6, 1.9 and 
2.4 pM). The results are sigmoidal curves as is shown in Fig. 1. It should be noted that the 
cytosolic mean Ca 2+ invariably increases with the level of [Ca 2+ ] a „ e . The vertical grey lines 
both in Fig. 1 and Fig. 2 are indicating the IP3 levels for which we have also studied several 
statistical properties of the puffs. Further insight can be obtained computing the effect of the 
IP3 over the mean fraction of open channels. The mean fraction of channel that is activated 
increase monotonically as the IP3 level increases following Hill curves (r 2 > 0.995 in all cases 
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as shown in the Fig. 2. The Hill curve is defined by 



where k is the half maximal value, V max is the maximum fraction of channels open, and is Hill 
coefficient. In the top panel of Fig. 3 we depict the Hill coefficient rih and in the bottom panel 
we depict k and V max for different values of [Ca 2+ ] aTO =. The Hill coefficient increases as [Ca 2+ ] a „ e 
increases. At physiological values of [Ca 2+ ] a „ e , the rih found in our simulation range from 1.15 
to 2.2 in very well agreement with the experimental values which range from 1.0 (Volpe et al., 
1990) to greater than 3.7 (Delisle, 1991). In the bottom panel of Fig. 3 we depict the others 
two fitted parameters V max and k. The maximum fraction of open channels V ma x ar e between 
[0.15,0.30] with a maximum at 0.80 of [Ca 2+ ] are , this range is slightly higher than the values 
found by Watras et al. (1991), around 0.14 in cerebellar cells (with = 1.3). It should be 
noted that the half-maximal k decrease with the level of [Ca 2+ ] ai)e and range from 0.10 to 0.7 
fiM IP3. The half-maximal found in the Watras experiment was 0.15 fiM IP3. 

Bezprozvanny et al. (1991) shown, in experiment where Ca 2+ channels reconstituted in 
lipid bilayers, that the fraction of open cannel as a function of log[Ca 2+ ] has an asymmetry 
bell-shaped curve, in experiment where calcium was clamped. The mathematical analysis of the 
master equation @ , allows to determine the open channel probability (to find the receptor in the 
state Xi^o) for a given [Ca 2+ ] and [IP3]. The steady-state distribution vector V st correspond 
to the normalized eigenvector associated to the eigenvalue zero of the equation of (j3J). Some 
algebraic manipulations lead to 



p si rrip 1 \c^w - KlK2 [IP3] [Ca ' J (*) 

n,i,o u^3j , i) i + ^ + R ^ [Ca2+] + K ^ R3 [i?3] [Ca2+]2 , , 



where K\ = k\/k^i, K2 = ^2/^-2 and K% = k^/k-^. Fig. 4 depicts the open channel probability 
for three different values of IP3 concentration as a function of [Ca 2+ ]. The result @ predicts the 
bell-shaped curve obtained experimentally by Bezprozvanny et al. (1991), with the symmetry 
axe correctly shifted to left. It should be noted the difference between Eq. and (jHJ). The 
first one can be understood as a open channel probability in physiological conditions described 
by Eq. (1-2). Thus the fraction of open channels is computed as the average of fraction of 
open channels over a period of time (10,000 seconds in this paper), where calcium is released 
(removed) to (from) cytosol many times. On the other hand, Eq. Q is a steady state probability 
for a fix level of calcium and IP3. This situation emulate the Bezprozvanny et al. experimental 
conditions, where the calcium was clamped. 

3.2 Statistical properties of puffs 

We now study the statistical properties of puffs obtained with the stochastic version of the OTM, 
the deterministic version of OTM is not able to address these properties. It is known that IP3 
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stimuli evokes repetitive spikes in the intracellular [Ca 2+ ] rather than simply increase [Ca 2+ ]. In 
the Fig. 5, we can see the typical spiking behavior of calcium release by a cluster with 20 channels 
(TV = 20) and some statistical properties of the observed puffs 10,000 seconds of simulations. 
In this case, we performed the simulation setting [Ca 2+ ] a „ e =1.6 /jM and IP3=0.35 /jM. At the 
top panel of Fig. 5 we depict the first 100 seconds of the temporal course of the intracellular 
calcium concentration. At the bottom-left panel, the scatter plot of IPIj versus IPF+i shown 
that there is not temporal correlation between two consecutive IPIs. At the bottom-center panel: 
scatter plot of puff amplitude versus puff duration shown that there is correlation between puff 
amplitudes and durations (r 2 =0.84). Small puff amplitudes are correlated with brief puffs, the 
correlation decreases as amplitude puff increases. Experimental observations of Thomas et al. 
(1998) indicate similar correlations, but with a smaller linear regression coefficient (r 2 =0.69). 
At the bottom-right panel we shown the puff size distribution, it should be noted that there are 
Ca 2+ releases of all sizes. All plots in the bottom panel of Fig. 5, and hereafter, were computed 
with puffs which amplitudes exceed the threshold 0.20 fiM., in order to avoid the background 
fluctuation. In this paper, the puff size is defined as the area between the [Ca 2+ ] time course 
and the threshold. 

It was experimentally established that when the IP3 concentration increases, the mean fre- 
quency of spiking raises, but the amplitude remain essentially constant. Fig. 6 displays the 
mean frequency of puffs as a function of the IP3 concentration for the some values of [Ca 2+ ] a „ e 
showed in Figs. 1-2 (0.8, 1.2, 1.6, 1.9 and 2.4 /jM.) obtained by stochastic simulations. We 
can observe that the onset of spiking behavior in response to the IP3 stimuli depends on the 
[Ca 2+ ] a „ e . For high levels of calcium, even low levels of IP3 are able to produce repetitive puffs. 
These responses are well characterized by Hill curves. The Hill coefficients associated to these 
curves were essentially independent of [Ca 2+ ] aw . level. For [Ca 2+ ] ai)e =2.4 fiM the Hill coefficient 
was 2.1 and decrease to 1.8 for [Ca 2+ ] a „ e =0.8 /xM. 

Furthermore, the number of channel recruited in the puffs varies stochastically from puff to 
puff. Important, and experimentally recorded, characteristics of these variabilities are amplitude, 
inter-puff interval, duration and size distributions (Bootman et al., 1997; Sun et al., 1998; 
Thomas et al., 1998; Callamaras and Parker 2000; Marchant and Parker 2001). The shape 
of the puff amplitude distribution depends on the concentration of IP3 and the [Ca + ] a t>e- 

In 

Fig. 7 we shown a diagram of the puff amplitude distributions in the [Ca 2+ ] a „ e -IP3 plane. In 
order to comparison the scale of the axis are the same for all plots in the diagram, between 
[0.20,1.80] fiM for the amplitudes of puffs (horizontal axe), and between [0,6000] for the number 
of events (vertical axe). No events (NE) were recorded for [IP3]= 0.07 fj,M and low level of 
[Ca 2+ ] a „ e =0.8 fjM). For lower [Ca 2+ ] a „ e , the amplitudes of the spontaneous puffs are typically 
smaller than 0.20 /jM and are regarded as fluctuations. For [Ca 2+ ] ai , e =0.8 fiM, or [IPs] = 0.07 
/xM, monotonically decreasing amplitude distribution are found. In the other cases a two-peak 
amplitude distributions are mainly found. The characteristic amplitudes depend on [Ca 2+ ] at)e , 
but not on [IP3]. For [Ca 2+ ] a „ e =1.6 /jM, the typical amplitude was around 0.80 /jM. The shape 
of the puff amplitude distributions shown in the Fig. 7 differ from those observed in experiments, 
which have one-peak shape (Sun et al., 1998; Thomas et al., 1998; Marchant and Parker, 2001). 
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Another important characteristic of the calcium dynamics is the distribution of time interval 
between two consecutive puffs. Recent experimental investigation with Xenopus oocytes has 
revealed that IPI distribution exhibits a single peak. In Fig. 8 we display another diagram 
in the [Ca 2+ ] a „ e -IP3 plane, but now with IPI distributions. In this diagram the range of the 
vertical axis are the same for all distributions, between [0,6000] events. While the ranges of 
IPI (horizontal axis) are between [0,80] seconds for the first column ([Ca 2+ ] ai)e =0.8 fjM.), and 
between [0,20] seconds for all other cases. For [Ca 2+ ] a „ e =0.8 /jM, the IPI distribution seems 
to decay exponentially. However, when the level of [Ca 2+ ] a „ e is increased, skewed single-peak 
distribution characterize the IPI distributions. The position of the single peak range from 1 to 
2 seconds increasing with the [Ca 2+ ] ai)e . The shape of the distributions shown in the Fig. 8 are 
consistent with observed by Marchant et al. (1999), but with a typical IPI smaller than that 
observed in the experiment (around 2.5 seconds). It is apparent that for high values of [IP3] and 
[Ca 2+ ] a „ e , the IPI are quite regular. 

We also want to illustrate the typical temporal behavior and some statistical properties, as 
we shown in the Fig. 5, for two different situations: i) high [IP3], low [Ca^J^g (Fig. 9); and 
ii) low [IP3], high [Ca 2+ ] a „ e (Fig. 10). It should be noted that different scale axis were used. 
In Fig. 9, where [IP3]=21 /jM, and [Ca 2+ ] a „ e =0.8 /jM., we found the distribution of puffs size 
decay exponentially, and there are very few great calcium release events (bottom-right panel). 
More broad distributions of IPI are yielded for this case than those shown in bottom-left panel 
of Figs. 5 and 10. There are also less IPI around 1-2 seconds (bottom-left panel) than those 
shown in Fig. 5 or in Fig. 10. The correlation between puff amplitude and duration seems to 
be similar than that shown in bottom-right panel of Fig. 5. When increase the [Ca 2+ ] a „ e until 
2.4 [iM and decrease the [IP3] to 0.07 fiM (Fig. 10), we found a slower exponential decay in 
the distribution of puffs size than in Fig. 9 (bottom-right panel). In this case the distributions 
of IPI (bottom-center panel of Fig. 10) looks like those shown in bottom-left panel of Figs. 5. 
Puff with small amplitude seems to be correlated with short duration. 



4 Discussion and conclusion 

We have developed a stochastic version of the Othmer-Tang model to study dynamical and 
statistical properties of Ca 2+ release of clusters of IP3-sensitive receptors. In comparison to 
others stochastic models (Falcke et al., 2000; Bar et al., 2000; Swillens et al., 1999), our receptor 
model is simpler and is represented by only four states. In this stochastic clustered IP3R model, 
the Ca 2+ diffusion between the cluster and the environment is ignored so that an isolated cluster 
can be discussed. However, the channels are assumed to be close enough and the instantaneous 
Ca 2+ diffusion within a cluster is so fast that the calcium concentration within a cluster is always 
homogeneous. 

We found that the model reproduce experimental result by Watras et al. (1991) where it is 
shown that the channel response to the IP3 follow a Hill curve. In intact cells, the long latency 
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and subsequent sigmoidal increase in cytosolic [Ca 2+ ] after addition of agonists (Berridge, 1997), 
are consistent with either positively cooperative activation of IP3RS by IP3 alone or by feedback 
from cytosolic [Ca 2+ ]. Distinguishing between these mechanisms, however, is difficult in intact 
cells, in which IP3 receptor activation inevitably increases cytosolic [Ca 2+ ]. In broken cells, the 
use of fluorescent indicators to measure rapid rates of Ca 2+ efflux is prone to similar difficulties. 
Results from a variety of tissues using superfusion (Finch et al., 1991), bilayer recording (Watras 
et al., 1991; Besprozvanny et al., 1991), fluorescent indicators (Champeil et al., 1989) or flash 
photolysis of caged IP3 in intact cells (Parker et al., 1996) suggest that responses to IP3 are 
either not positively cooperative (Finch et al., 1991; Parker et al., 1996; Watras et al., 1991) 
or only slightly so (Besprozvanny et al., 1991; Champeil et al., 1989). Each approach has its 
limitations and together the results provide no clear indication of whether responses to IP3 are 
positively cooperative in the absence of positive feedback from cytosolic Ca 2+ . In our numerical 
simulation we found that, even in absence of cooperativity in the model, the Hill coefficient 
increase from 1.15 to 2.2 when [Ca 2+ ] ai)e increase. However, the model is not able to explain 
higher values of as that reported by Marchant and Taylor (1997). This results suggest that 
feedback from cytosolic Ca 2+ plays a key role in the channel response to IP3, but to explain 
high values of found in some tissues, a cooperative mechanism could be mandatory. 

Other important aspect is the frequency encoding. It is well known that IP3 stimuli evoke 
repetitive spikes in the intracellular Ca 2+ concentration rather than simply elevating the level 
of Ca 2+ . Moreover, when the [IP3] is increased, the mean frequency of spiking raises, but the 
amplitude remains essentially constant. Thus a continuous level of stimuli signal is converted 
into a frequency-encoded signal: the number of spikes per time unit (Tang and Othmer, 1995). 
A common problem found in this kind of converter devices is the competition between two 
desirable goals: i) high sensitivity, the system ideally should be able to detect even weak signals; 
and ii) large dynamic range, the system should not saturate over various orders of magnitude 
of input intensity. For this reason the Hill coefficient (when the response is well-fitted by a Hill 
curve) its desirable to be less than 1. In the results derived from our numerical simulation, we 
can observe that the mean frequency of the spiking behavior in response to the IP3 stimuli is well 
characterized by Hill curves. The Hill coefficient, between 1.8-2.1, was essentially independent 
of [Ca 2+ ] a „ e . This value is too high given narrow dynamical range in order to frequency-encode 
the IP3 signal. However, this narrow dynamical range observed in isolated cluster, could be 
translates into a wider dynamical range as result of the collective phenomenon by considering 
an intercluster communication mechanism. This interesting discussion is the subject of our 
current research and will be discussed in a forthcoming paper. 

The small number of the IP3RS in a cluster introduce stochastic oscillation into the calcium 
release dynamics. Experimental data suggest that the localized calcium release varies in a con- 
tinuous fashion due to stochastic variation in both numbers of channels recruited and durations 
of channel openings. It is not obvious which aspects of Ca 2+ puffs are originally determined by 
the dynamics of the Ca 2+ channels, which properties are determined by the diffusion and Ca 2+ 
binding kinetics, and which properties are induced from the measurement procedures (Cheng et 
al., 1999; Izu et al., 1998). Although in several experiments one-peak amplitude distributions 
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were observed, Cheng et al., (1999) suggested that the original calcium puffs should have an 
exponential decaying decreasing amplitude distribution. Our stochastic simulations shown that 
the amplitude distributions varies in a continuous fashion from exponential decay to two-peak 
distributions. Even for a small threshold (as 0.20 //M) a peak is observed in the amplitude 
distributions for high values both of [IP3] and [Ca 2+ ] ai)e . However, the puffs size (i.e. ampli- 
tude x duration) distributions have essentially an exponential decay distribution over the range 
studied here. These results indicate that a fixed puff morphology which is sometimes assumed 
in literature (Pratusevich and Balke, 1996; Izu et al., 1998; Cheng et al., 1999) is not a good 
assumption for Ca 2+ puff analysis. Furthermore, in our markovian simulation, we found that 
the amplitude of puff is correlated with its duration in agreement with experimental observation 
reported by Thomas et al., (1998). 

Other useful way to characterize the calcium puffs dynamics is using the IPI distributions. 
The stochastic simulation of the OTM shown that the IPI distributions varies from exponential 
decay to a skewed bell-shape distributions. For physiological values of IP3 and [Ca 2+ ] ai)e the 
distributions have a maximum around 1.5 seconds. The shape is consistent with distribution 
observed by Marchant et al. (1999), however, the typical IPI observed in the experiment, 
around 2.5 seconds, was slightly higher. In contrast to the deterministic models, we shown 
that two consecutive IPI are not correlated. For this reason, the random spiking behavior or 
stochastic oscillations have an essentially different nature that the periodic oscillations generated 
by deterministic models. 

We want to remark that the stochastic oscillations are different from the stochastic excitabil- 
ity discussed by Keizer and Smith (1998). For the stochastic excitability, once [Ca 2+ ] randomly 
becomes larger than a threshold, a fast release (action-potential-like) of [Ca 2+ ] followed by a 
refractory period can be observed. For the stochastic oscillation studied here there is not such 
a threshold. More broad distributions of puff amplitudes and IPI are yielded for stochastic 
oscillation. 
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Tables and Figures 



Kinetic parameters 
name value 


Other parameters 
name value 


fci 12.0 QuMx s)- 1 

k-i 8.0 s- 1 

k 2 23.4 (fiMx s)- 1 

k- 2 1-65 s- 1 

k 3 2.81 (/iMx s)" 1 

A:_ 3 0.21 s- 1 


v r 1.85 

70 0.1 s- 1 

71 20.5 s- 1 

pi 8.5 /iMx s _1 
p 2 0.065 /M 
iV 20 



Table 1: Parameters values for the OTM. 
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gure 1: Mean value of [Ca 2+ ] versus IP3 concentration, for six different concentrations of the 
erage intracellular calcium [Ca 2+ ] ai , e , computed over 10,000 seconds. 
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Figure 2: Fraction of open channels versus IP3 concentration, for six different concentrations of 
the average intracellular calcium [Ca 2+ ] me , computed over 10,000 seconds. 
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Figure 3: Top panel: Hill coefficient for IP3-induced calcium release versus [Ca 2+ ] me , obtained 
fitting data displayed in figure 2. Bottom panel: k and V max values of the Hill curves for IPs- 
induced calcium release versus [Ca 2+ ] ai)e , obtained fitting data displayed in Fig. 2. 
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Figure 4: The theoretically predicted probability to find the open channel as a function of 
cytosolic [Ca 2+ ] for three different values of IP3 concentration (0.01 fiM. solid line, 0.1 /iM. 
dotted line, and 1/dVI dot-dashed line). 
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Figure 5: Top panel: 100 seconds of the stochastic temporal behavior of the calcium release from 
a cluster with 20 IP 3 Rs (N = 20, [Ca 2+ ] at)e =1.6 and IP 3 =0.35). Bottom-left panel: scatter plot 
of IPI(z) versus IPI(z + l) shown that there is not temporal correlation between two consecutives 
puffs. Bottom-center panel: scatter plot of puff amplitude versus puff duration shown that there 
is some correlation between puff durations and amplitudes of puffs. Bottom-right panel: puff 
area distribution shown that there is Ca2+ release of all sizes. All plots in the bottom panel were 
computed with puffs which amplitudes exceeds the threshold 0.20 //M, recorded when simulate 
10,000 seconds. 
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Figure 6: Mean frequency versus IP3 concentration, for five different values of [Ca 2+ ]ave- All 
plots were computed with puffs which amplitudes exceeds the threshold 0.20 /xM, recorded when 
simulate 10,000 seconds. 
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Figure 7: Diagram of the puff amplitude distributions is shown in the [Ca 2+ ] at , e -IP3 plane. The 
distributions was computed with all puffs, which amplitudes exceeds the threshold 0.20 /dVt, 
recorded in f 0,000 seconds. In this condition no events (NE) were recorded for [IPa] = 0.07 
/xM and [Ca 2+ ] a „ e = 0.8 jjM. The range of the axis are the same for all distributions, between 
[0.20,1.8] for the amplitudes of puffs (horizontal axe), and between [0,6000] for the number 
of events (vertical axe). 
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Figure 8: Diagram of the inter puff interval (IPI) distributions is shown in the [Ca ^] a ^e" IPs 
plane. The distributions was computed with all puffs, which amplitudes exceeds the threshold 
0.20 fiM, recorded in 10000 seconds. In this condition no events (NE) were recorded for [IPs] = 
0.07 fjM and [Ca 2+ ] a „ e = 0.8 fM. The range of the vertical axis are the same for all distributions, 
between [0,6000] events. While the ranges of IPI (horizontal axis) are between [0,80] seconds for 
the first column ([Ca 2+ ] at , e =0.8 /xM), and between [0,20] seconds in all other cases. 
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Figure 9: Top panel: 100 seconds of the stochastic temporal behavior of the calcium release 
from a cluster with 20 IP 3 Rs (N = 20, [Ca 2+ ] aw; =0.8 and IP 3 =21.0 fM). Bottom-left 
panel: scatter plot of IPI(i) versus IPI(i + 1) shown that there is not temporal correlation 
between two consecutives puffs. Bottom-center panel: scatter plot of puff amplitude versus puff 
duration shown that there is some correlation between puff durations and amplitudes of puffs. 
Bottom-right panel: puff area distribution shown that there is Ca 2+ release of all sizes. All plots 
in the bottom panel were computed with puffs which amplitudes exceeds the threshold 0.20 //M, 
recorded in 10,000 seconds. 
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Figure 10: Top panel: 100 seconds of the stochastic temporal behavior of the calcium release 
from a cluster with 20 IP 3 Rs (N = 20, [Ca 2 +] ave =2A fiM and IP 3 =0.07 fM). Bottom-left panel: 
scatter plot of IPI(i) versus IPI(i + 1) shown that there is not temporal correlation between two 
consecutives puffs. Bottom-center panel: scatter plot of puff amplitude versus puff duration 
shown that there is some correlation between puff durations and amplitudes of puffs. Bottom- 
right panel: puff area distribution shown that there is Ca 2+ release of all sizes. All plots in 
the bottom panel were computed with puffs which amplitudes exceeds the threshold 0.20 /jlM, 
recorded in 10,000 seconds. 
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